Simulations of frustrated Ising Hamiltonians using quantum approximate optimization

Novel magnetic materials are important for future technological advances. Theoretical and numerical calculations of ground-state properties are essential in understanding these materials, however, computational complexity limits conventional methods for studying these states. Here we investigate an alternative approach to preparing materials ground states using the quantum approximate optimization algorithm (QAOA) on near-term quantum computers. We study classical Ising spin models on unit cells of square, Shastry-Sutherland and triangular lattices, with varying field amplitudes and couplings in the material Hamiltonian. We find relationships between the theoretical QAOA success probability and the structure of the ground state, indicating that only a modest number of measurements (≲100) are needed to find the ground state of our nine-spin Hamiltonians, even for parameters leading to frustrated magnetism. We further demonstrate the approach in calculations on a trapped-ion quantum computer and succeed in recovering each ground state of the Shastry-Sutherland unit cell with probabilities close to ideal theoretical values. The results demonstrate the viability of QAOA for materials ground state preparation in the frustrated Ising limit, giving important first steps towards larger sizes and more complex Hamiltonians where quantum computational advantage may prove essential in developing a systematic understanding of novel materials. This article is part of the theme issue ‘Quantum annealing and computation: challenges and perspectives’.


Introduction
Quantum magnetism has been a major focus in condensed matter research, driven by the potential for new disruptive applications ranging from quantum computing to quantum sensing [1]. Quantum material properties are intrinsically related to the structure of the ground states. However, exact ground states are notoriously challenging to calculate classically, requiring the field to resort to using semi-classical limits [2][3][4][5] or fully quantum approaches with restricted applicability [6][7][8][9][10][11][12]. New, fully quantum computational tools are required to understand current problems including frustrated two-dimensional quantum magnets currently explored by bulk neutron scattering and thin film susceptibility [13,14]. Digital and analogue quantum simulators have emerged as a new tool for the simulation of quantum many-body phenomena towards efficient modelling of exotic quantum phases of matter beyond classical tractability [15,16]. They are naturally suited for magnetic Hamiltonians since spins can be directly mapped to qubits. Non-trivial phases in magnetic systems, such as frustrated phases [17], spin glasses [18], and topologically ordered phases [19,20], have been realized on multiple qubit platforms using a variety of techniques.
In this paper, we investigate an alternative approach to preparing materials ground states using the quantum approximate optimization algorithm (QAOA) [21] on near-term quantum computers. We apply QAOA to lattices of interest in materials science, considering the classical Ising limit (equivalently, S = ∞) where standard QAOA is directly applicable. This serves as a stepping stone towards truly quantum problems such as the XY and Heisenberg models in the fully frustrated limit, which will require further algorithmic research and modifications to the approach presented here. Our results validate that QAOA achieves sufficient accuracy for the simpler classical limit and provides insights into algorithmic behaviour for material lattice problems.
We consider lattice instances with varying degrees of frustration. The smallest building block of a frustrated magnetic Hamiltonian is an anti-ferromagnetic triangular motif of three spins where all the bonds cannot be satisfied simultaneously. In these materials, exchange interactions compete such that it is impossible to satisfy them all simultaneously (figure 1). If all spin configurations are equally favourable, frustration can lead to non-ordered states such as spin liquids [5], spin glasses [22] or plaquette states [8], each with distinct signatures.
We solve three different types of Hamiltonians for unit cells pictured in figure 2. The first is a square unit cell Hamiltonian, which exhibits only simple ferromagnetic and antiferromagnetic phases in the infinite size (thermodynamic) limit. The second is the celebrated Shastry-Sutherland lattice. Interestingly, this problem already lends itself to materials applications and experimental data analysis. Among other examples, it is conjectured to describe the class of rare-earth tetraborides (ErB 4 , TmB 4 and NdB 4 ) and allows a direct comparison with several existing results both theoretical [23][24][25][26][27] and experimental [28][29][30][31][32][33]. The third case is the more complex Ising triangular lattice which represents a maximally frustrated problem with an infinite number of possible ground states in the infinite size limit [34,35]. We compute theoretical probabilities to prepare the ground state for each of these 9-spin Hamiltonians under varying choices of the external field and coupling parameters and compare these theoretical results against computations on a trapped ion quantum computer.  We choose N = 9 spins as this is the logical minimum number of spins required to construct a unit cell of the Shastry-Sutherland lattice and also it is a feasible size for the Quantinuum quantum computer (with 12 qubits available at the time this work was completed). We focus on p = 1 layers of the QAOA algorithm; for larger N instances more layers p of QAOA will be needed to maintain a significant success probability [36][37][38]. Implementations on quantum computers will also have to overcome predicted limitations due to noise [39][40][41][42] including an exponential scaling in the number of measurements with circuit size, depth, and other factors [43]. Noise in modern quantum computers has negative consequences for all quantum algorithms, not only QAOA. Ongoing testing and development of these devices is necessary to assess realistic performance scaling in the presence of noise and to determine whether QAOA, or other quantum algorithms, will ultimately succeed in providing a useful computational advantage over conventional approaches.

Ising Hamiltonian and model unit cell lattices
A single unpaired spin on the outermost orbital of a magnetic ion constitutes a s = 1/2 state which is implemented straightforwardly on a qubit. In a magnetic material, several such spins in a lattice interact via pairwise superexchange interactions J α . The nature and strength of these interactions, J α , depend on several factors. These include the distance between the magnetic ions (typically J α scales as the inverse cube of the distance between the ions), the shape of the orbitals, the symmetry of the lattice and the local crystal fields. Magnetic frustration can arise in the lattice, for example, if spins arranged on a triangular motif in the lattice experience equal Ising anti-ferromagnetic interactions. Such a magnetic frustration can arise via a combination of straight edge bonds J 1 which are either horizontal or vertical, and the diagonal bonds J 2 . This is given by the Hamiltonian where the first sum is over the nearest neighbours (NN), the second sum is over the diagonal next-nearest neighbours (NNN), and s = (s 1 , . . . , s N ) lists the spin orientations s i ∈ {1, −1} of the N spins. We study anti-ferromagnetic couplings with positive J 1 , J 2 . The term h represents a longitudinal magnetic field (parallel to the spin axis), which for the real material represents either a mean crystal field or an external magnetic field. The unit cell motif of the Hamiltonian is shown in figure 2. Materials described by this model are being actively researched in condensed matter physics. The Ising Shastry-Sutherland model is a special case of a model, inspired by the geometry of real materials, where some but not all of the diagonal bonds are present (figure 2b).
The triangular lattice is shown in figure 2c. In all cases, we consider open boundary conditions. Analytical ground state properties of Ising models on Shastry-Sutherland and triangular lattices have been derived analytically [44,45]. Multiple methods have been proposed to solve for ground states of Ising Hamiltonians and related optimization problems, notably Integer Programming method [46], Simulated Annealing [47] and its variants, Large Neighbourhood Search [48] and Quantum or Quantum-inspired physical annealing devices. Among them, the Integer Programming method solves exactly but suffers from exponential scaling of computational time. Simulated Annealing and Large Neighbourhood Search are both herustics methods that promise faster runtime but there is no guarantee of the solution qualities. Quantum annealers, digital annealers and coherent Ising machines are hardwares dedicated to solving Ising models [49,50]. Depending on the connectivity of these various (qu)bits, every backend could be good at a different task-parallel tempering machines could be good at finding the classical phases [26], while others could reveal intricate dynamical behaviour in a transverse field Ising universality [51]. For frustrated lattice problems, QAOA allows us to sample different ground states with certain probabilities due to quantum randomness, whereas classical and deterministic algorithms may generate a ground state efficiently, but fail to explore the states which might arise because of a coherent superposition between all the spins. Additionally, future extensions beyond the Ising limit also become an exciting possibility.

(a) Ground state magnetization phase diagrams
We consider the nine-spin unit cells with geometries in figure 2a-c which represent the number of spins required to simply construct a unit cell of the Shastry-Sutherland lattice. In materials represented by Bravais lattices, these unit cells repeat periodically to realize the very large lattices in a real material. We computed ground states for each unit cell by evaluating (2.1) for each possible spin configuration to identify the lowest energy states, for varying choices of h and J 2 , with J 1 = 1 taken as the unit of energy. We plot the magnetization  Starting with the non-frustrated square lattice with J 1 -only interactions with simple ferromagnetic or antiferromagnetic ground states, the degree of frustration is tuned progressively by (i) addition of J 2 bonds and (ii) bringing J 2 → J 1 . The triangular lattice with uniform coupling parameters represents the maximally frustrated limit with highly degenerate solutions. The Shastry-Sutherland model represents a scenario with the minimum number of J 2 bonds required to realize a fully frustrated lattice. The solution of these states represents a problem of polynomial time complexity in two dimensions and without a magnetic field.
The ground state for a given h shows a number of magnetization plateaus, where each plateau has a different proportion of spins pointing up. Unsurprisingly, at large h, the ground state for each lattice is ferromagnetic, in regions C, G and E in figure 2d-f, respectively. The situation becomes more interesting at small h in regions 'A', the ground state is anti-ferromagnetic with magnetization M = 1/9, as five spins are aligned with the field while four are anti-aligned. For fields 8/3 ≤ h ≤ 4 and small J 2 , there is a ground state with M = 7/9 in which a single spin in the centre of the unit cell is anti-aligned with the field. Besides frustration, these states are also determined by the finite size of the unit cells, where the central spin is distinguished as the only spin with four interactions in the square lattice. As J 2 and h are varied, frustration leads to a variety of different ground states for the Shastry-Sutherland and triangular lattices, with varying magnetizations in figure 2e,f, with ground states in electronic supplementary material. These are true ground states of the 9-spin Hamiltonians, with boundary spins playing a big role. In the infinite size limit, we expect the ground states to be progressively less dependent on the boundary, and more on the symmetry of the Hamiltonian, which we discuss in the next section.

(b) Finite size effects
The finite sizes of our lattice unit cells, as well as the unusual M = 7/9 ground state noted in the previous section, raise questions of how the ground states for our unit cells match with ground states that would be obtained in the large size limit, and the minimum number of spins that are needed to achieve quantitative behaviour consistent with large sizes. To address these questions, we computed the magnetization of triangular and Shastry-Sutherland lattices of N = n × n spins to analyse the size-dependent behaviour. Owing to the exponential complexity of the problem, we used neal [52], a software implementation of simulated annealing to approximate the ground state configurations with h, J 2 ∈ [0, 6] and J 1 = 1. Each combination of h and J 2 was run 50 times and the solution with minimum Ising energy was picked. Examples with 3 ≤ n ≤ 30 are shown in figure 3. We assessed convergence of the global phase diagrams to the large size limit by computing the root mean square error (RMSE) between the target lattice's and 30×30 lattice's magnetization across points in the phase diagram. We fitted the RMSE to both a power-law as well as an exponential form (see electronic supplementary material, figure S1), and we find that the power-law scaling with n exponent γ and prefactor a fits the RMSE better. The equation takes the form:    result based on an overall RMSE demonstrates that a 15 × 15 spin grid is already obtaining results close to the much larger 30×30 grid. Based on this trend, we expect that finite size effects in M will diminish quickly with the size of the lattice, indicating that lattices of only a few hundred spins may diminish the errors sufficiently to achieve a realistic 'bulk', and therefore meaningful results for comparison with experiments which probe bulk properties, such as diffraction and heat-capacity. This suggests that quantum processors with hundreds of qubits, achievable within the noisy intermediate-scale quantum era [55], may be capable of meaningful applications for materials science applications.

Quantum approximate optimization algorithm
Quantum computers offer a route to overcoming issues associated with identifying ground states through conventional methods. One approach to address these problems uses the quantum approximate optimization algorithm, which was originally designed to find approximate solutions to difficult combinatorial optimization problems [21] that are often expressed in terms of Ising Hamiltonians [56]. Empirical performance of QAOA has been characterized for a variety of combinatorial problems [36,[57][58][59][60] and this has also led to generalizations [61][62][63][64][65][66] that have been applied to preparing chemical ground states [67] as well as ground state preparation for onedimensional [38,68] and two-dimensional [69] quantum spin models in theory and experiment [70].
To formulate our Ising problems in a structure that is suitable for QAOA, we first express the Ising Hamiltonian (2.1) in terms of a quantum Hamiltonian operator where H(z) comes from (2.1) taking s i = 1 − 2z i for each component |z i in the total basis state |z . This gives an encoding of the Ising spin problem (2.1) that is useful for QAOA, where we will sample eigenstates |z to try to identify the ground state of the Ising problem. To find solutions, QAOA uses a quantum state prepared with p layers of unitary evolution, where each layer alternates between Hamiltonian evolution under the Ising Hamiltonian H and under a 'mixing' where the initial state |ψ 0 = 2 −N/2 z |z is the ground state of −B represented in the computational basis. The parameters γ = (γ 1 , . . . , γ p ) and β = (β 1 , . . . , β p ) are typically chosen to minimize the expectation value of the energy H , though other objectives have also been studied [67,71,72]. The minimization is typically accomplished using a quantum-classical hybrid feedback loop, shown schematically in figure 5. For a given set of parameters γ and β, a set of states |ψ p (γ , β) is prepared and measured by a quantum computer. The measurement results are sent to a conventional (classical) computer to compute the classical objective function. If the objective function is not converged relative to previous evaluations, then the conventional computer uses an optimization routine to select new parameters γ , β . The process is repeated until convergence to a minimal objective with optimized parameters γ * , β * . The final result is taken as the measurement result |z * that gives the lowest energy H(z * ). In the best case, z * = z ground is a ground state, while more generally z * may be a low-energy state that is an approximate solution to the problem. An analytic proof has demonstrated that QAOA can prepare an exact ground state |z ground of Ising Hamiltonians as p → ∞ [21,56]. Apart from the formal proof of convergence at large p, there has been significant interest in applying QAOA at small p, where approximations exceeding conventional lower bounds have been observed in simulations [36,73] and predicted for large problems in certain contexts [74]. Realizing such advantages on devices with hundreds of qubits or more is an important topic of ongoing research as quantum computing technologies continue to develop.
For the materials problems we consider, we are interested in preparing the ground states of the Ising Hamiltonian. We compute exact ground states for our unit cells in figure 2 by evaluating all eigenvalues of the Hamiltonian using equations (2.1) and (3.1) to identify the lowest energy state. For some phases the number of ground states is N ground > 1 due to degeneracies, while for other phases there is a single ground state N ground = 1, as pictured in electronic supplementary material. To assess QAOA performance, we compute the average ground-state probability where the sum contains a single term in the case of a non-degenerate ground state or multiple terms in the degenerate case. Analytically, the probabilities are given by the Born rule P(z) = | z|ψ p (γ , β) | 2 , while for experiments on a quantum computer they are given by the frequencies of measurement results, P(z) = N(z)/N tot , where N(z) is the number of times |z was measured and N tot is the total number of measurements. If QAOA identifies a ground state then |z * = |z ground and P ground > 0, while if QAOA only finds sub-optimal solutions then P ground = 0. Ground-state preparation is a goal specific to the materials problem context we are interested in here. This is, importantly, a departure from the standard goal of QAOA in the context of approximate combinatorial optimization, where the goal is to find approximate solutions that are not necessarily the ground states. While QAOA is not expected to efficiently find exact ground states for generic NP-hard optimization problems, it may still prove useful for finding ground states of specific structured problems such as materials problems on a lattice similar to those we explore here [38,[68][69][70]. We use numerical calculations to assess the theoretical performance of QAOA for ground-state preparation. These demonstrate the ideal performance of QAOA in exact pure state calculations that use matrix multiplication to evaluate (3.3). This gives an ideal baseline for later comparison against results from a noisy quantum computer, where errors lead to mixed states with degraded performance. We use p = 1 QAOA layers throughout this section and our results.
To identify QAOA states for our Ising problems, we must determine optimized QAOA parameters γ * 1 and β * 1 . We choose regions to evaluate parameters in determining γ * 1 and β * 1 as follows. QAOA is periodic when β 1 → β 1 ± π [57], hence we consider −π/2 ≤ β 1 ≤ π/2, which gives all unique β 1 up to symmetries. The periodicity of the γ 1 parameter is more complicated, as it depends on the Hamiltonian in exp(−iγ 1 H) in (3.3). Here we focus on γ 1 intervals near the origin and dependent on the magnitude of the Hamiltonian terms, which has been highly successful in previous work [75]. The basic idea is that the QAOA unitary exp(−iγ 1 H) changes at varying speeds with respect to γ 1 , depending on the Hamiltonian coefficients J 1 , h and J 2 . When the Hamiltonian coefficients increase, then γ 1 should decrease to obtain a similar unitary. The rate at which the unitary changes with respect to γ 1 is related to the average magnitude of the Hamiltonian coefficients where E NN is the number of nearest-neighbour interactions and E NNN is the number of nextnearest-neighbour interactions. Previous work on generic Ising Hamiltonians with h = 0 has shown that high-quality solutions are obtained at small γ 1 with an empirical scaling of optimized parameters similar to γ * 1 ∼ 1/ι. The scaling 1/ι compensates for the varying rates of evolution that are present for varying choices of the Hamiltonian, and also limits the interval of γ 1 values that are explored, simplifying the optimization [75]. Based on these ideas, we choose γ 1 in the interval −0.55 × π/ι ≤ γ 1 ≤ 0.55 × π/ι.
We  We have found in searches over much larger γ 1 intervals that the local optima for P ground and H do not always approximately align as in figure 6, which can lead to poor P ground at optimized H in these larger intervals. However, this does not appear to be a prevalent issue for the smaller ι-dependent γ 1 intervals in cases we have looked at. The results are somewhat sensitive to the specific choice of γ 1 interval, however, our choice of −0.55 × π/ι ≤ γ 1 ≤ 0.55 × π/ι gives satisfactory results across the varying lattices.
To identify optimized parameters, we perform a grid search over γ 1 and β 1 for each Hamiltonian considered. We evaluate the QAOA states in (3.3) on 201 evenly spaced intervals with −π/2 ≤ β 1 ≤ π/2 and over 300 evenly spaced intervals with −0.55 × π/ι ≤ γ 1 ≤ 0.55 × π/ι for a total of 201 × 300 = 60 300 grid evaluations for each Hamiltonian. This approach gives optimal parameters in our intervals up to coarse graining in the grid search. We select parameters γ * 1 and β * 1 that optimize H . The optimized parameters γ * 1 and β * 1 do not necessarily give the optimal ground state probabilities that are possible from QAOA. The reason is that the average energy optimization accounts for energies and probabilities of all states, which together may yield low energies at parameter choices that are not optimal for the ground states alone [36]. To assess performance, we further compare our parameter choices against parameters γ * 1 and β * 1 that directly optimize P ground . The direct optimization of P ground is used here for benchmarking purposes and is not a realistic approach for large problems where the ground states are unknown. For our small problems, the comparison gives an idea of how the ground state probabilities from a standard optimization of H compare against the best ground state probabilities that could be obtained by QAOA in our set-up.

(b) Quantum computations of QAOA
We next investigate the performance of QAOA using the Quantinuum H1-2 quantum computer. H1-2 contains trapped-ion qubits and uses lasers to implement gates on these qubits. Typical error rates are reported as 3.5 × 10 −3 for two-qubit gates and 1 × 10 −4 for single-qubit rotation gates [76]. In addition to the device H1-2, we also use the H1-1E device 'emulator' to simulate noisy device behaviour. This gives results that approximately correspond to expected device behaviour while avoiding the financial expense and wait times that are associated with running the device. The emulator models a variety of device-specific noise processes for the H1-class computers, including depolarizing noise, leakage errors, crosstalk, dephasing in transport and qubit idling errors [77].
We test QAOA on the H1-2 using the QCOR software stack [78]. The QCOR stack translates the series of unitary operators expressing QAOA into quantum circuits for H1-2; see electronic supplementary material, appendix B, for details. The QCOR program used for submitting jobs to the device as well as our calculations are available online, cf. 'Data accessibility'.
Furthermore, modern quantum computers are known to be affected by state preparation and measurement (SPAM) errors as well as gate infidelities from a variety of physical sources. We assessed SPAM errors expected in our quantum computations using the device emulator, with details in electronic supplementary material, appendix C. The probability to observe no error in circuits we tested was approximately 96%, with errors distributed approximately uniformly across qubits. We account for these errors using an independent bit-flip model and associated SPAM matrixP, which transforms an ideal set of measurement results to the expected noisy set of results. The inverse matrixP −1 can then be applied to our noisy measurements from the quantum computer to approximately correct for SPAM errors. A technical issue arises in that the mitigated measurement probabilities can sometimes be negative, due to the approximate nature of the mitigation scheme. This leads to a second mitigation scheme that additionally sets all negative probabilities to zero and renormalizes so the total probability is one. We use each of these approaches to attempt to correct the small SPAM errors we expect from the quantum computer, as described in detail in electronic supplementary material.

Results
In this section, we consider the results from QAOA applied to the materials lattices of figure 2. We take J 1 = 1 as the unit of energy and analyse the success of QAOA in preparing ground states at variable h and J 2 , first in numerical simulations ( §3a) and then in quantum computations on a trapped-ion quantum computer ( §3b).

(a) Ground-state measurement probabilities
We first consider theoretical probabilities to measure the ground state with QAOA and how these vary for different parameter choices in the Hamiltonian. We begin with the simple square lattice in figure 2a,d, which does not exhibit frustration as there are no triangles in the interaction graph. The probability to measure the ground state for varying h is shown in figure 7. Figure 7a shows the probabilities obtained from optimizing the standard objective H , while figure 7b shows the bestcase results based on optimizing P ground , as described in §3a. The probabilities in each case are similar, demonstrating that optimizing H is nearly as successful in increasing the ground-state probability as a direct optimization. The average ground-state probability shows distinct behaviours for each of the three ground states at varying h, visually separated by dotted lines. In the anti-ferromagnetic ground state at small h, the probability P ground approximately oscillates between h = 0 and h = 2, with small probabilities observed near integer values of h and larger probabilities near h = 1/2 and h = 3/2. The M = 7/9 ground state with 8/3 < h < 4 has a near-constant probability of ≈ 0.06. At h = 4 the ground state becomes ferromagnetic and P ground increases significantly, with monotonic increases at larger h. We rationalize the varying success probabilities in the figure as attributable to structures of the ground states at varying h and the interplay with the structure of the QAOA state in (3.3). We show in electronic supplementary material, appendix D that QAOA can exactly prepare the ferromagnetic ground state when h J 1 for arbitrary lattice sizes, based on the fact that exp(−iγ H) ≈ exp(−iγ h N i=1 Z i ) in this limit. This is consistent with the behaviour in the figure, where P ground increases monotonically with h for the ferromagnetic ground state at h ≥ 4. We further show in electronic supplementary material, appendix E how the anti-ferromagnetic ground state probability is maximized at h = 0.5, and we devise large γ 1 parameters that can further improve these results. (We did not include larger γ 1 parameters in our numerical searches as this can lead to poor P ground at parameters that optimize H for the frustrated lattices, as  remarked in §3a.) However, the mechanism for anti-ferromagnetic ground state preparation here depends on the specific choice of the 3 × 3 lattice, and it is not clear how QAOA will behave for other lattice sizes. For the M=7/9 phase, the QAOA state is more complicated, as it is a superposition of many basis states that depends on the optimized parameters, and we do not have an analytic account for this behaviour. The optimized parameters that create each QAOA state in the square-lattice phase diagram are shown in electronic supplementary material, appendix F.
Ground-state probabilities P ground for the Shastry-Sutherland and triangular lattices are pictured in figures 8 and 9, respectively. Ground-state probabilities from optimizing H are presented in panels (a) while panels (b) show the best case probabilities from a direct optimization of P ground as described in §3a. These probabilities show patterned behaviour, with distinct probabilities P ground observed throughout most of each individual region A,B,. . ., with significant differences in P ground between different regions. At small J 2 , there are oscillations in the probability for preparing the anti-ferromagnetic ground states at small h, and large success probabilities for the ferromagnetic ground state at large h, as foreshadowed by results from the square lattice. On the other hand, as the J 2 coupling increases, the triangular and Shastry-Sutherland lattices experience increased frustration, with competing interactions within the triangular motifs in figure 2. The average ground-state probability decreases significantly as J 2 increases and frustration becomes dominant. This is especially evident when J 2 h, for example in the top left of each of figures 8b and 9b. The P ground are mostly uniform across h and J 2 within each region, qualitatively similar to the nearly uniform probability for the M = 7/9 state at varying h for the square lattice in figure 7. Ground-state probabilities are typically 0.01, indicating that only 100 measurements are expected to identify a ground state. We now turn to computations on a trapped-ion quantum computer, to benchmark and assess performance of QAOA on a real quantum computing device.

(b) QAOA quantum computations
Here we assess QAOA performance in preparing ground states on a trapped-ion quantum computer. Ultimately, our aim is to validate the idea that a current quantum computing technology is capable of preparing each ground state of a frustrated Ising Hamiltonian using QAOA. An important first step is to assess whether optimized parameters from our theoretical calculations are also optimized for the device, or whether further optimization is needed to determine device-specific optimized parameters.    values and we attribute this to noise in the device emulator, which introduces errors that cause the energy to deviate from its ideal minimum value. Despite these errors, the shape of the landscape is similar to our theoretical calculations, with best performance observed near the theoretical parameters that minimize H , and energies that tend to increase away from these parameters.
We further validate that the H1-2 trapped-ion device itself is consistent with the emulator in the data points that begin with 'H1-2'. These actually yield better energies than the emulator, and are within one standard error of the mean from our theoretical results. The results from the device and emulator indicate that the energy landscape as a function of the QAOA parameters γ 1 and β 1 is consistent between our theoretical calculations, the quantum device, and emulator. We therefore proceed with our theoretically optimized parameters to evaluate success in ground-state preparation using the quantum computer.

(ii) Quantum calculations of ground states
We now perform calculations on the Honeywell H1-2 quantum computer to analyse success probability in ground-state preparation. We consider points in each region of the Shastry-Sutherland lattice, using parameters listed in table 1 that correspond to local minima in H , similar to the minima used to evaluate theoretical performance in §4a). We post-process the measurement results using the SPAM mitigation model with probabilities P(z) ≥ 0 (see §3b), to give a minor correction to the observed results that is designed to counteract this known source of error. Figure 11 shows  For a closer comparison of probabilities, we plot error bars denoting the theoretical standard error of the mean s.e.m.= P(z)(1 − P(z))/N shots . The s.e.m. defines a range in which we expect about two-thirds of estimated P est (z) = N(z)/N shots are expected to be found, where N(z) is the number of measurement results of a given ground state and N shots is the total number of measurements in table 1. The probabilities from the quantum computation are largely consistent with the pure state results, with the majority of results within one s.e.m. from the ideal P(z), as expected in finite sampling to estimate the ground state probability. There are only two large deviations for states 'k' and 'i', which may be related to noise in the device. The quantum computations succeed in preparing ground states with probabilities comparable to pure-state expectations.

Conclusion
We analysed QAOA as an approach for preparing materials ground states on three types of Ising Hamiltonians with longitudinal magnetic fields, focusing on nine-spin unit cells as a starting size that is amenable to simulations and calculations on a trapped-ion quantum computer.
We applied QAOA to the nine-spin Ising unit cell problems to assess its success in ground state preparation. We found that the theoretical success probability depends significantly on the structure of the ground state, while it is mostly insensitive to the precise Hamiltonian parameters, which can vary within regions consistent with a fixed ground state. Each Hamiltonian yields a ferromagnetic ground state in the presence of large magnetic fields, and QAOA achieved large success probabilities for these relatively simple states. The probabilities for other types of ground states were more variable, and tended to decrease as next-nearest-neighbour couplings became stronger with associated frustration in the lattice. For all of these nine-spin states, we typically find success probabilities indicating that 100 measurements are expected to be necessary from an ideal quantum computer for these problem instances.
To assess QAOA performance under realistic conditions, we implemented the algorithm on a trapped-ion quantum computer. These quantum computations succeeded in observing each of the 17 ground states of the Shastry-Sutherland unit cell. The quantum computations yielded ground state probabilities that were consistent with theoretical expectations based on pure states, indicating that noise was not a significant issue at the sizes and depths tested. This suggests that calculations with current technology can likely be extended to greater QAOA depth parameters p, and to larger sizes as greater numbers of qubits become available. At greater depths and sizes we expect higher performance and more realistic results in comparison with the large size limit, respectively. While the ground states and associated phase diagrams for our nine-spin unit cells were found to have significant finite size effects relative to the large-size limit, the errors from finite size effects on the classically calculated magnetization phase diagrams on n × n lattices up to n = 30 was found to be suppressed rapidly with n, with small errors at n = 15 indicating that only hundreds of spins may be necessary to reproduce large scale behaviour. This provides a baseline of hundreds of qubits for quantum computational experiments that seek to explain materials science problems, which may be accessible to near-term quantum computers in coming years.
Assessing scaling of the ground-state probability with size N will be an essential aspect of extending this approach to larger sizes N. This includes numerical simulations to quantify how the ground-state probabilities depend on the number of spins N and number of QAOA layers p; previous works have shown p ∝ N maintains a large ground-state probability ( 0.7) for simple models in different contexts [37,38], but future work is needed to test scaling in the current model. Benchmarking on quantum computers is also essential to understand how real noise processes effect scalability.
Thus from our results we envision QAOA can be successfully applied to somewhat larger lattice problems as quantum computing technologies develop and larger number of qubits become available. These could be used for optimization of Ising lattice problems as we have here, with increasing sizes that potentially describe real bulk properties of materials in the N → ∞ limit. Additionally, the QAOA algorithm is general, with the application of a magnetic field, and hence could be explored for lattice problems which are NP-Hard [79]. But a more promising future direction leveraging the full benefit of this approach is to extend and modify QAOA to prepare ground states of quantum Hamiltonians such as the XY and Heisenberg models, which can lead to a variety of quantum phenomena not captured in the Ising model, such as quantum spin glasses [80], spin nematicity [81], Berzinski-Kosterlitz-Thouless states [82,83] and long-range entangled states such as Dirac string excitations [84], the likes of which exist in twodimensional frustrated quantum spin liquids and spin ice. Many of these topics are fiercely researched and are of considerable interest and importance for future quantum technologies and devices. Conventional numerical methods for understanding these states are hindered by the exponential size of the Hilbert space, making it difficult to generate a theoretical understanding of experimental observations. QAOA or related generalizations [38,61,[67][68][69][70] offer a potential route to overcome conventional computing bottlenecks. Some successes along these lines have been observed in certain contexts, however, advances in methodology and quantum computing technologies are needed to extend these methods to complicated and larger-scale problems where quantum computational approaches may have a significant impact in understanding and developing materials for technological applications.
The data are provided in electronic supplementary material [85].